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Abstract. A Bernoulli free boundary problem with geometrical constraints is studied. The domain £1 
is constrained to lie in the half space determined by x\ > and its boundary to contain a segment of the 
hyperplane {x\ = 0} where non-homogeneous Dirichlet conditions are imposed. We are then looking for 
the solution of a partial differential equation satisfying a Dirichlet and a Neumann boundary condition si- 
• multaneously on the free boundary. The existence and uniqueness of a solution have already been addressed 

and this paper is devoted first to the study of geometric and asymptotic properties of the solution and then 
| to the numerical treatment of the problem using a shape optimization formulation. The major difficulty and 

' originality of this paper lies in the treatment of the geometric constraints. 
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1 Introduction 

Let (0, x\, ...,xn) be a system of Cartesian coordinates in R with N > 2. We set = {M, N : x\ > 0}. 
Let K be a smooth, bounded and convex set such that K is included in the hyperplane {x\ = 0}. We define 
a set of admissible shapes O as 

O = {$7 open and convex, K C dQ}. 
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Figure 1 : The domain f2 in dimension two. 



We are looking for a domain f2 E 0, and for a function u : $7 
system 



R such that the following over-determined 



(1) 
(2) 
(3) 
(4) 



-An = in tt, 
u = 1 on K, 
u = ondn\K, 



|Vn|=l on r := (dQ \ K) n R 



N 



has a solution; see Figure[T]for a sketch of the geometry. Problem CD)-© is afree boundary problem in the 
sense that it admits a solution only for particular geometries of the domain fi. The set Y is the so-called free 
boundary we are looking for. Therefore, the problem is formulated as 



(5) 



(J 7 ) : Find VL G O such that problem (Q]) — © has a solution. 



This problem arises from various areas, for instance shape optimization, fluid dynamics, electrochemistry 
and electromagnetics, as explained in |[T] [8l [TUl [TTJ . For applications in N diffusion, we refer to ll26l and for 
the deformation plasticity see 0. 

For our purposes it is convenient to introduce the set L := (dQ \ K) n {x\ = 0}. Problems of the type 
(J 7 ) may or may not, in general, have solutions, but it was already proved in [24] that there exists a unique 
solution to (J 7 ) in the class O. Further we will denote Q* this solution. In addition, it is shown in [24] that 
dQ* is C 2+a for any < a < 1, that the free boundary 50* \ K meets the fixed boundary K tangentially 
and that L* = (dft* \K)C\ {x\ = 0} is not empty. 

In the literature, much attention has been devoted to the Bernoulli problem in the geometric configura- 
tion where the boundary d£l is composed of two connected components and such that $7 is connected but not 
simply connected, (for instance for a ring-shaped 0), or for a finite union of such domains; we refer to [3] [9] 
for a review of theoretical results and to (4J |T3l H31 HSJ [22j for a description of several numerical methods for 
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these problems. In this configuration, one distinguishes the interior Bernoulli problem where the additional 
boundary condition similar to (@]) is on the inner boundary, from the exterior Bernoulli problem where the 
additional boundary condition is on the outer boundary The problem studied in this paper can be seen as a 
"limit" problem of the exterior boundary problem described in [ 16], since <9S7 has one connected component 
and is simply connected. 

In comparison to the standard Bernoulli problems, (J 7 ) presents several additional distinctive features, 
both from the theoretical and numerical point of view. The difficulties here stem from the particular geomet- 
ric setting. Indeed, the constraint Q C is such that the hyperplane {xi = 0} behaves like an obstacle 
for the domain Q and the free boundary dfl \ K. It is clear from the results in ll24l that this constraint will 
be active as the optimal set L* = (dQ* \ K) n {xi = 0} is not empty. This type of constraint is difficult to 
deal with in shape optimization and there has been very few attempts, if any, at solving these problems. 

From the theoretical point of view, the difficulties are apparent in [24], but a proof technique used for 
the standard Bernoulli problem may be adapted to our particular setup. Indeed, a Beurling's technique and 
a Perron argument were used, in the same way as in lfT6l[l7l[T8l . 

Nevertheless, the proof of the existence and uniqueness of the free boundary is mainly theoretical and 
no numerical algorithm may be deduced to construct T. From the numerical point of view, several problems 
arise that will be discussed in the next sections. The main issue is that T is a free boundary but the set 
L = (dfl \ K) n {x\ = 0} is a "free" set as well, in the sense that its length is unknown and should be ob- 
tained through the optimization process. In other words, the interface between L and T has to be determined 
and this creates a major difficulty for the numerical resolution. 

The aim of this paper is twofold: on one hand we perform a detailed analysis of the geometrical prop- 
erties of the free boundary T and in particular we are interested in the dependence of T on K. On the other 
hand, we introduce an efficient algorithm in order to compute a numerical approximation of Q. In this way 
we perform a complete analysis of the problem. 

First of all, using standard techniques for free boundary problems, we prove symmetry and monotonic- 
ity properties of the free boundary. These results are used further to prove the main theoretical result of the 
section in Subsection 13.31 where the asymptotic behavior of the free boundary, as the length of the subset 
K of the boundary diverges, is exhibited. The proof is based on a judicious cut-out of the optimal domain 
and on estimates of the solution of the associated partial differential equation to derive the variational for- 
mulation driving the solution of the "limit problem". Secondly, we give a numerical algorithm for a 
numerical approximation of Q. To determine the free boundary we use a shape optimization approach as 
in lfl3l[T4"l[l9l . where a penalization of one of the boundary condition in ([T])-© using a shape functional is 
introduced. However, the original contribution of this paper regarding the numerical algorithm comes from 
the way how the "free" part L of the boundary is handled. Indeed, it has been proved in the theoretical study 
presented in ll24l that the set L = (dQ\K) n {xi = 0} has nonzero length. The only equation satisfied on 
L is the Dirichlet condition, and a singularity naturally appears in the solution at the interface between K 
and L during the optimization process, due to the jump in boundary conditions. This singularity is a major 
issue for numerical algorithms: the usual numerical approaches for standard Bernoulli free boundary prob- 
lems |[T3l [T4l [T9l l20l cannot be used and a specific methodology has to be developed. A solution proposed 
in this paper consists in introducing a partial differential equation with special Robin boundary conditions 
depending on an asymptotically small parameter e and approximating the solution of the free boundary 
problem. We then prove in Theorem [3] the convergence of the approximate solution to the solution of the 
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free boundary problem, as e goes to zero. In doing so we show the efficiency of a numerical algorithm that 
may be easily adapted to solve other problems where the free boundary meets a fixed boundary as well as 
free boundary problems with geometrical constraints or jumps in boundary conditions. Our implementation 
is based on a standard parameterization of the boundary using splines. Numerical results show the efficiency 
of the approach. 

The paper is organized as follows. Section [2] is devoted to recalling basic concepts of shape sensitivity 
analysis. In Section [Sj we provide qualitative properties of the free boundary T, precisely we exhibit sym- 
metry and a monotonicity property with respect to the length of the set K as well as asymptotic properties 
of r. In Section HI the shape optimization approach for the resolution of the free boundary problem and a 
penalization of the p.d.e. to handle the jump in boundary conditions are introduced. In section [51 the shape 
derivative of the functionals are computed and used in the numerical simulations of (J 7 ) in Section [6] and |7J 

2 Shape sensitivity analysis 

To solve the free boundary problem (J 7 ), we formulate it as a shape optimization problem, i.e. as the mini- 
mization of a functional which depends on the geometry of the domains f2 C O. In this way we may study 
the sensitivity with respect to perturbations of the shape and use it in a numerical algorithm. The shape 
sensitivity analysis is also useful to study the dependence of on the length of K, and in particular to 
derive the monotonicity of the domain Q* with respect to the length of K. 

The major difficulty in dealing with sets of shapes is that they do not have a vector space structure. In 
order to be able to define shape derivatives and study the sensitivity of shape functionals, we need to con- 
struct such a structure for the shape spaces. In the literature, this is done by considering perturbations of an 
initial domain; see I6l fl5ll27l . 

Therefore, essentially two types of domain perturbations are considered in general. The first one is a 
method of perturbation of the identity operator, the second one, the velocity or speed method is based on 
the deformation obtained by the flow of a velocity field. The speed method is more general than the method 
of perturbation of the identity operator, and the equivalence between deformations obtained by a family of 
transformations and deformations obtained by the flow of velocity field may be shown (6l|27|. The method 
of perturbation of the identity operator is a particular kind of domain transformation, and in this paper the 
main results will be given using a simplified speed method, but we point out that using one or the other is 
rather a matter of preference as several classical textbooks and authors rely on the method of perturbation 
of the identity operator as well. 

For the presentation of the speed method, we mainly rely on the presentations in (6j [271 • We also 
restrict ourselves to shape perturbations by autonomous vector fields, i.e. time-independent vector fields. 
Let V : ~R N — > II be an autonomous vector field. Assume that 

(6) V G V k (R N ,R N ) = {V G C k (R N , R N ), V has compact support}, 

with k > 0. 

For t > 0, we introduce a family of transformations T t (V)(X) = x(t, X) as the solution to the ordinary 
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differential equation 

(7) f J t x ^ x ) = V(x(t,X)), 0<t<r, 

\ x(0,X) = X£R N . 

For r sufficiently small, the system (O has a unique solution 11271 . The mapping T t allows to define a family 
of domains Q t = T t (V)(Q) which may be used for the differentiation of the shape functional. We refer to 
(U Chapter 7] and (27 , Theorem 2.16] for Theorems establishing the regularity of transformations T t . 

It is assumed that the shape functional J($7) is well-defined for any measurable set Q, C IR^. We 
introduce the following notions of differentiability with respect to the shape 

Definition 1 (Eulerian semiderivative). Let V G V k (R N ,M, N ) with k > 0, the Eulerian semiderivative of 
the shape functional J(£l) at £1 in the direction V is defined as the limit 

, 8 ) ww-^mzm, 

when the limit exists and is finite. 

Definition 2 (Shape Differentiability). The functional J(Q) is shape differentiable (or differentiable for 
simplicity) at if it has a Eulerian semiderivative at Q in all directions V and the map 

(9) V^dJ(n,V) 

is linear and continuous from ^(R^, R ) into R The map © is then sometimes denoted X7J(Q) and 
referred to as the shape gradient of J and we have 

(10) dJ(tl,V) = (VJ(il),V r ) :D -fc( R JV jR JV )jI ,fc( R JV jR JV ) 

When the data is smooth enough, i.e. when the boundary of the domain Q and the velocity field V 
are smooth enough (this will be specified later on), the shape derivative has a particular structure: it is 
concentrated on the boundary <9S1 and depends only on the normal component of the velocity field V on the 
boundary dQ,. This result, often called structure theorem or Hadamard Formula, is fundamental in shape 
optimization and will be observed in Theorem 0] 



3 Geometric properties and asymptotic behaviour 

In shape optimization, once the existence and maybe uniqueness of an optimal domain have been obtained, 
an explicit representation of the domain, using a parameterization for instance usually cannot be achieved, 
except in some particular cases, for instance if the optimal domain has a simple shape such as a ball, ellipse 
or a regular polygon. On the other hand, it is usually possible to determine important geometric properties 
of the optimum, such as symmetry, connectivity, convexity for instance. In this section we show first of all 
that the optimal domain is symmetric with respect to the perpendicular bissector of the segment K, using a 
symmetrization argument. Then, we are interested in the asymptotic behaviour of the solution as the length 
of K goes to infinity. We are able to show that the optimal domain f2* is monotonically increasing for the 
inclusion when the length of K increases, and that converges, in a sense that will be given in Theorem 
12 to the infinite strip (0,1) x R. 

The proofs presented in Subsections 13. II and 13.21 are quite standard and similar ideas of proofs may be 
found e.g. in flBKHHTTKIll. 
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3.1 Symmetry 

In this subsection, we derive a symmetry property of the free boundary. The interest of such a remark is 
intrinsic and appears useful from a numerical point of view too, for instance to test the efficiency of the 
chosen algorithm. 

In the two-dimensional case, we have the following result of symmetry: 

Proposition 1. Let Q* be the solution of the free boundary problem (Q])-©. Assume, without loss of gener- 
ality that (Oxi) is the perpendicular bissector of K. Then, f2* is symmetric with respect to (Ox\). 

Proof. Like often, this proof is based on a symmetrization argument. It may be noticed that, according to 
the result stated in [24, Theorem 1], Q is the unique solution of the overdetermined optimization problem 



(So): 

where 



minimize J(0, u) 

subject to fl G O, u G H(Q), 



fl-(fi) = {u£ H 1 (fi),u = 1 on -u = on dQ\Kand\Vu\ = lonT}, 

and 

J(tt,u) = / \Vu(x)\ 2 dx. 



From now on, Q* will denote the unique solution of (Bo), K being fixed. We denote by SI the Steiner 
symmetrization of Q with respect to the hyperplane X2 = 0, i.e. 

SI = jcc = (x',x 2 ) such that - ^|0(a/)| < x 2 < ^\Q(x')\,x' G ft'j , 



where 



and 



Q' = {x' G R such that there exists x 2 with (x 1 , x 2 ) G fl*} 



U(x') = {x 2 G R such that (x',x 2 ) G Q}, x' G W. 
By construction, is symmetric with respect to the (Ox\) axis. Let us also introduce u, defined by 

u : x S H- sup{c such that x G oj*(c)}, 

where co*(c) = {x G Q* : u(x) > c}. Then, one may verify that u G H(Q) and Polya's inequality (see 
ifBl) yields 

j(h,u) < j(n*,u*). 

Since (fl*, u*) is a minimizer of J and using the uniqueness of the solution of (Bo), we get fl* = Q. □ 

Remark 1. This proof yields in addition that the direction of the normal vector at the intersection ofT and 
(Oxi) is (Ox\). 



6 



3.2 Monotonicity 



In this subsection, we show that fl* is monotonically increasing for the inclusion when the length of K 
increases. For a given a > 0, define K a = {0} x [—a, a}. Let {F a ) denote problem {F) with K a instead of 
K and denote fl a and u a the corresponding solutions. We have the following result on the monotonicity of 
fl a with respect to a. 

Theorem 1. Let < a < b, then fl a C ^lb- 
Proof. According to (24), (F a ) has a solution for every a > and <9r2 a is C 2+a , < a < 1. We argue by 
contradiction, assuming that fl a <f_ fib- Introduce, for t > 1, the set 

fit = {x £ fl a : tx G ^a}- 

We also denote by iv~ t := {0} x [—ta,ta] and T t := 9f2t\(cK7j n (OX2)). The domain fl t is obviously a 
convex set included in fl a for £ > 1. Now denote 

t min := inf{t > 1, fl t C fi;,}. 

On one hand, Sl a C is equivalent to i m j n = 1. On the other hand, if fl a <f_ fib, then t m i n > 1 and 
for t large enough, we clearly have fl t C fib, therefore t m i n is finite. In addition, if fi a ^ fib we have 
r trai „ n r 6 / 0. Now, choose y € r tmin n r 6 . Let us introduce 

U trm. n -X eU t \-^ U a (t min x). 

Then, ut verifies 

—Aut . =0 in fi+ , 
iu, ■ = 1 on iv"+ , 
ut =0 on T+ , 

so that, in view of fit min C fib and K tmin C -fTb, the maximum principle yields > ut min in fit min . 
Consequently, the function h = Ub — ut min is harmonic in fit m< „, and since h(y) = 0, h reaches its lower 
bound at y. Applying Hopf's lemma (see 0) thus yields d n h{y) < so that \Vub(y)\ > |Vu tmin (y)|. 
Hence, 

1 = \Vu b (y)\ > \Vu tmin (y)\ = t min > 1, 
which is absurd. Therefore we necessarily have t m i n = 1 and fl a C fib. □ 



3.3 Asymptotic behaviour 

We may now use the symmetry property of the free boundary to obtain the asymptotic properties of fl a when 
the length of K goes to infinity, i.e. we are interested in the behaviour of the free boundary T a as a — > 00. 

Let us say one word on our motivations for studying such a problem. First, this problem can be seen as 
a limit problem of the "unbounded case" studied in lfl8l Section 5] relative to the one phase free boundary 
problem for the p-Laplacian with non-constant Bernoulli boundary condition. Second, let us notice that the 
change of variable x' = x/a and y' = y/a transforms the free boundary ©-((U) problem into 

(11) - Az 

(12) z 

(13) z 

(14) |V*| 



= in fl, 

= 1 on K\ , 

= ondfl\K h 

= a onT = (<9fi \ #1) n R.; 
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which proves that the solution of (fTTT)-([T4l) is ^i/ o (0 a ), where h x i a denotes the homothety centered at the 
origin, with ratio I /a. Hence such a study permits also to study the role of the Lagrange multiplier associated 
with the volume constraint of the problem 

J mmC(Q) where C(O) = min f n \Vun\ 2 , u = 1 on K%, u = on d£l\K\} 
\ VL quasi-open, |f2| = m, 

since, as enlightened in lfT31 Chapter 6], the optimal domain is the solution of (fTTV(fl4l for a certain con- 
stant a > 0. The study presented in this section permits to link the Lagrange multiplier to the constant m 
appearing in the volume constraint and to get some information on the limit case a — > +oo. 

We actually show that T a converges, in an appropriate sense, to the line parallel to K a and passing 
through the point (1,0). Let us introduce the infinite open strip 

S =]0, l[xE, 

and the open, bounded rectangle 

R(b) =]0,l[x] -6,6[C S. 

Let 

us ■ x G S i— > 1 — x\. 

Observe that, since Q a is solution of the free boundary problem ([T])-©, the curve r a n { — cl ^ X2 ^ o} is 
the graph of a concave C 2,a function X2 h-> ^(^2) on [—a, a]. We have the following result 

Theorem 2. The domain Vt a converges to the strip S in the sense that for all b > 0, we have 

(15) ip a — > 1, uniformly in [—b, b], as a — ^ +00. 
We also have the convergence 

u a — > us in -ff 1 (i?(6)) as a — > 00, 

for the solution u a of CL])-©. 
Proof. Let us introduce the function 

v a {xi) = u a (x 1 ,0). 

According to [15, Proposition 5.4.12], we have for a domain Q of class C 2 and u : £1 — > R of class C 2 

(16) An = A v u + Ud n u + d 2 u, 

where Apu denotes the Laplace-Beltrami operator. Applying formula (fT6l ) in the domains 

w (c) := {x £ n a ,u a (x) > c}, 
we get Arn a = on du a (c), Au a = due to (Q} and thus 

(17) d 2 n u a = -U a d n u a on duj a (c), 

where n is the outer unit normal vector to ui a {c) and Ti a (x) denotes here the curvature of du a (c) at a point 
x G du a (c). Thanks to the symmetry of Q, a with respect to the xi-axis, we have d n u a (x\, 0) = v' a (x\) and 
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d%u a (xi,0) = v"(x\) for xi > 0. According to [24], the sets uj a (c) are convex. Therefore ~H a is positive 
on dco a (c) and v a (x\) is non-increasing. Thus 

(18) vZ(x 1 ) = -H a (x 1 ,0)v' a (x 1 )>0, 

which means that v a is convex. Let m a be such that T a n (R x {0}) = (m a , 0), i.e. the first coordinate of 
the intersection of the xi-axis and the free boundary T a . The function v a satisfies 

(19) -<(xi)<0 far a* G]0,m a [, 

(20) v a (0) = 1, 

(21) Da(m a ) = 0, 

(22) <(m a ) = -1. 

In view of (fT9l >. t> a is convex on [0, m a }. Since t> a (0) = 1 and v a (m a ) = 0, then 

v a {xi) < 1 - — • 

Furthermore, m a < 1, otherwise, due to the convexity of v a , the Neumann condition (1221 would not be 
satisfied. Since S7 a is convex, this proves that f2 a C S and that Q a is bounded. 

Moreover, from Theorem [T] the map a \— > S7 a is nondecreasing with respect to the inclusion. It follows that 
the sequence (m a ) is nondecreasing and bounded since Q a C S. Hence, (m a ) converges to < 1. 
Let us define 

Uoo(xi) = 1 • 

moo 

The previous remarks ensure that for every a > 0, v a < Uoq. 

Let T>(a) be the line containing the points (0,a) and (t/) a (b),b) and T(a) the line tangent to T a at 
(i^a(b), b). Let sx>(a) and sr(a) denote the slopes of T>(a) and T(a), respectively. For a fixed 6 G (0, a), 
we have 

6 - a 

SDia) = , — >■ — oo as a — )• oo, 

since < ^ a < 1. Due to the convexity of O a , we also have sj-{a) < sx>(a). Therefore 

sj-(a) — > — oo as a —7- oo. 

Thus, the slopes of the tangents to F a go to infinity in Q a n R(b). Furthermore, due to the concavity of the 
function tp a , we get, by construction of T>{a), 

-^(a - x 2 ) < ip a (x2) < rnoo, Va > 0, Vx 2 G [-6, £>]. 
Hence, we obtain the pointwise convergence result: 

(23) lim 4> a (x2) = rrioo, Vx 2 G [-b,b], 

a— >+oo 

which proves the uniform convergence of ip a to moo as a — )• +oo. 

From now on, with a slight misuse of notation, u a will also denote its extension by zero to all of S. 
Finally, let us prove the convergence 



in iT 1 (i? 0O (6)), as a — > oo, 
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where Roo(b) denotes the rectangle whose edges are: Si = {0} x [—6,6], S2 = [0,moo] x {6}, S3 = 
{moo} x [-6,6] and S 4 = [0,moo] x {-6}. 

According to the zero Dirichlet conditions on S3 and using Poincare's inequality, proving the H 1 -convergence 
is equivalent to show that 

(24) / |V(« a - ^oo)| 2 ->■ asa->oo. 

For our purposes, we introduce the curve £2(0) described by the points X a ^ solutions of the following 
ordinary differential equation 

(25) { ^(*) = Vu a (X a>b (t)), t > 0, 

I X a>b (0) = (0,b). 

The curve £2(0) is naturally extended along its tangent outside of Q a . £2(0) can be seen as the curve 
originating at the point (0, 6) and perpendicular to the level set curves of Q a . We also introduce the curve 
£4(0), symmetric to £2(0) with respect to the xi-axis. £4(0) is obviously the set of points Y a ^ solutions 
of the following ordinary differential equation 

(26) J ^(t) = V« a (^W), t>0, 

1 F o>6 (0) = (0,-6). 

Then the set Q a (b) is defined as the region delimited by the X2-axis on the left, the line parallel to the X2-axis 
and passing through the point (moo, 0) on the right and the curves £2(0) and S 4 (a) at the top and bottom. 
We also introduce the set S3 (a) := Q a (b) D ({moo} x R). See Figure |2]for a description of the sets i?oo(6) 
and Q a (b). 

Since Roo(b) C Q a (b) (see Figure [2]), we have 

/ \V{ Ua - Uoo )\ 2 < I |V(u a -Uoo)| 2 . 

jRa^b) JQa(b) 

Using Green's formula, we get 

/ |V(u a -Moo)| 2 = / |V(u a - Uoo)\ 2 + / IV^a-Uoo)! 2 

jQa{b) JQ a (b)nn a JQa(b)\Qa 



(U a - Uoo)A(u a - Uoo) - (u a - Uoo)A(u a - Uoo) 

i(b)nn a JQ a (b)\n a 



+ / (u a - Uoo)d u (u a - Uoo) \ 

JdQaib) + JTo 



(U a - Uoo)d n ±(u a - Uoo), 

ldQ a (b) " Jr a nQ a (b) 

where v denotes the outer normal vector to Q a {b) on the boundary dQ a {b), n is the outer normal vector 
to Q, a on the boundary T a and d n ± is the normal derivative on T a in the exterior or interior direction, the 
positive sign denoting the exterior direction to Q a . The functions u a and Uoo are harmonic, and using the 
various boundary conditions for u a and Uoo we get 

/ \V(u a - Uoo)\ 2 = I _ (u a - Uoo)d u (u a - Uoo) + / Uoo- 

JQ a (b) iE 2 (fl)US 4 (a) JT a nQ a (b) 
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Figure 2: The sets Roo(b) and Q a (b). 



According to (T23T ) and using Uqo = on S3, we get 



/ u OQ = u 00 (ifj a (x2))y / l + ^' a (x 2 ) 2 dx2 -> as a ->• 00, 

Jr a nQ a (fe) 

where we have also used the fact that ij)' a {x2) - > for all a?2 € [—6, 6]. The limit function depends only 
on x\, thus we have dyU^ = on E2 U £4. Denote now ip a : [0, moo] — > K- the graph of £2(0) (which 
implies that —ip a is the graph of £4(0)). The slope of the tangents to the level sets of u a converge to —00 
as a — > 00 in a similar way as for T a , therefore d X2 u a {x\, ij) a {xi)) converges uniformly to in [0, moo] as 
a — > 00, and in view of d25l > we have that ip a — >■ b uniformly in [0, moo] and since (u a — u^o) is uniformly 
bounded in Q a we have 

(27) / _ (u a - u^dyUoo as a -> 00. 

is 2 (a)US 4 (a) 

In view of the definition of Q a {b), the outer normal vector v to Q a (b) at a given point on £2(0) U £4(0) 
is colinear with the tangent vector to the level set curve of J7 a passing though the same point. Therefore 
d v u a = on £2(0) U £4(0) and we obtain finally 

(28) 0< / \V{u a - Uoo )\ 2 < I |V(u a -Uoo)| 2 ^0 as a ^00. 

The end of the proof consists in proving that moo = !• Let us introduce the test function ip as the 



11 



solution of the partial differential equation 

r -Av9 = in Q a (b)^ 
(29) < tp = onSiUS 2 (a)US 4 (a) 

[ cp = l on £ 3 (a). 

It can be noticed that 99 E i/ 1 (i2 QO (6)). 

Using Green's formula and the same notations as previously, we get 



/ V(u a - Uoo) • Vlf = / V(li a - Uoo) • V99 + / V (u a - Uoo) ■ V ip 

JQ a (b) JQ a (b)nn a JQ a (b)\n a 

/ <fA(u a - Uoo) ~ / ¥^{ u a - U 

jQ„.(b)r\n a JQ„.(b)\n a 



^00 



<pd v {u a - Uoo) + I ipd n (u a - Uoo) - / V 3 

E 2 (a)UE 4 (a) JS 3 ((i) Jr a nQ a (fe) 



'S 3 (a) m oo 7r Q nQ a (6) 
According to (|23T ). and since we deduce from (l28l) that 



/ V(n a — Uqo) • V99 — ^ as a — > 00, 

JQa{b) 



we get 



which leads to 

1 



— " i <P = o, 

E 3 (a) m oo Js 3 (a) 

1 ] |E 3 (a)| =0. 



moo 

In other words, m^o = 1, which ends the proof. □ 



4 A penalization approach 
4.1 Shape optimization problems 

From now on we will assume that N = 2, i.e. we solve the problem in the plane. The problem for N > 2 
may be treated with the same technique, but the numerical implementation becomes tedious. A classical 
approach to solve the free boundary problem is to penalize one of the boundary conditions in the over- 
determined system £[])-© within a shape optimization approach to find the free boundary. For instance one 
may consider the well-posed problem 

(30) -Aui = in O, 

(31) ui = l onK, 

(32) ui = on dn \ K. 
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and enforce the second boundary condition (01) by solving the problem 

/ yj \ f minimize J(f2) 

(33) {Bl): \ subject to n G0, 

with the functional J defined by 

(34) j(n) = J (d nUl + i) 2 dr. 

Indeed, using the maximum principle, one sees immediately that u\ > in SI and since u\ = on dQ. \ K, 
we obtain d n u\ < on <9S1 \ K. Thus |Vui| = —d n u\ on dQ, \ K and the additional boundary condition 
dU) is equivalent to d n u\ = — 1 on Y. Hence, (1341) corresponds to a penalization of condition (@]). On one 
hand, if we denote u\ the unique solution of (Q])-((4]) associated to the optimal set Q*, we have 

J(Q*) = 0, 

so that the minimization problem (l33l has a solution. On the other hand, if J(Q,*) = 0, then | Vn^| = 1 on 
T and therefore it^ is solution of CQ)-©- Thus (J 7 ) and (£>i) are equivalent. 

Another possibility is to penalize boundary condition ([3]> instead of (0]) as in (B\), in which case we 
consider the problem 

(35) -A« 2 = in 0, 

(36) u>2 = 1 on if, 

(37) U2 = on L, 

(38) <9 n u 2 = -l onT, 



minimize J(Q) 
subject to SI G 0, 



and the shape optimization problem is 

(39) (B 2 ) : 
with the functional J defined by 

(40) J(O) = / {u 2 fdT. 



Although the two approaches (B\) and (B2) are completely satisfying from a theoretical point of view, it 
is numerically easier to minimize a domain integral rather than a boundary integral as in (l34l and (l40l . 
Therefore, a third classical approach is to solve 

.... ,„ . f minimize J(S1) 

(41) I subject to flVo, 

with the functional J defined by 

(42) j(n)= / K- U2 ) 2 . 
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Figure 3: Polar coordinates with origin Aj, and such that 0j = corresponds to the semi-axis tangent to F. 

For the standard Bernoulli problems (3]|9l, solving (B3) is an excellent approach as demonstrated in lfT3l 
[14j[19l. However, we are still not quite satisfied with it in our case. Indeed, it is well-known that due to the 
jump in boundary conditions at the interface between L and F in (|37T)-(|38T). the solution U2 has a singular 
behaviour in the neighbourhood of this interface. To be more precise, let us define the points 

{a 1 ,a 2 } :=Xnr, 

and the polar coordinates (r,, 6i) with origin the points Aj, i = 1,2, and such that Oi = corresponds to 
the semi-axis tangent to F; see Figure [3] for an illustration. Then, in the neighbourhood of A4, 112 has a 
singularity of the type 

s i {n,e i ) = c{A i )^F i cos(e i /2), 

where c(Ai) is the so-called stress intensity factor (see e.g. Ifl2~ll2ll ). 

These singularities are problematic for two reasons. The first difficulty is numerical: these singularities 
may produce inacurracies when computing the solution near the points {Ai, A2}, unless the proper numer- 
ical setting is used. It also possibly produces non-smooth deformations of the shape, which might create in 
turn undesired angles in the shape during the optimization procedure. The second difficulty is theoretical: 
since F is a free boundary with the constraint Q C R+ , the points {Ai, A2} are also "free points", i.e. their 
optimal position is unknown in the same way as F is unknown. This means that the sensitivity with respect 
to those points has to be studied, which is doable but tedious, although interesting. The main ingredient in 
the computation of the shape sensitivity with respect to these points is the evaluation of the stress intensity 
factors c(Ai). 



4.2 Penalization of the partial differential equation 

In order to deal with the aforementionned issue, we introduce a fourth approach, based on the penalization 
of the jump in the boundary conditions (l37l)-(P38l) for u%. Let e > be a small real parameter, and let 
ip e G C(Et + ,R + ) be a decreasing penalization function such that ip £ > 0, tjj £ has compact support [0,/3 e ], 
and with the properties 

(43) p e -> as e -> 0, 

(44) Ve(0) -> ooase -> 0, 

(45) ^eOl) -> Oase -> 0, Vxi > 0. 

A simple example of such function is given by 



(46) M^i) =e _1 (max(l -e _9 xi,0)) 2 l 



E+) 
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with q > 0. Note that ip £ is decreasing, has compact support and verifies assumptions (|43T)-(|4"51). with 
(3 £ = e q . We will see in Proposition [2] that the choice of tp £ is conditioned by the shape of the domain. Then 
we consider the problem with Robin boundary conditions 

(47) -Au 2 ,e = in 0, 

(48) u 2 , e = 1 on if, 

(49) d n u 2 ^ + i>e{x 1 )u 2 ,e = -1 on 5ft \ if. 

The function u 2j£ is a penalization of u 2 in the sense that u 2>£ — > u 2 as e — > in H 1 ^) if ip £ is properly 
chosen. The following Proposition ensures the H 1 -convergence of u 2:£ to the desired function. It may be 
noticed that an explicit choice of function ip £ providing the convergence is given in the statement of this 
Proposition. 



Proposition 2. Let ft be an open bounded domain. Then for ip £ given by (1461 ), there exists a unique solution 
to (l47l-d49l) which satisfies 

(50) u 2>£ — > u 2 in H l (Q) as e — >■ 0. 

Proof. In the sequel, c will denote a generic positive constant which may change its value throughout the 
proof and does not depend on the parameter e. 
We shall prove that the difference 

V £ = U 2 - U 2)E . 

converges to zero in H 1 ^!). The remainder v £ satisfies, according to (|35T)-(|38T) and (l47"T)-(|4"9l) 

(51) -Av e = in ft, 

(52) v £ = on K, 

(53) d n v £ + ip £ {0)v £ = 1 + d n u 2 on L, 

(54) d n v £ + ip £ (xi)v £ = *p £ (xi)u 2 on T. 

Multiplying by v £ on both sides of (IBTT ). integrating on ft and using Green's formula, we end up with 

(55) / \Vv £ \ 2 + [ (v £ ) 2 ^ £ = [ ^ £ u 2 v £ + [ (l + d n u 2 )v £ . 
Jn Jan Jr Jl 

Since v £ = on K we may apply Poincare's Theorem and (|55T ) implies 

(56) ^Klll^n) ^ c (ll^ u 2||L2 (r) ||v e || L 2 (r) + \\l + d n u 2 \\ L 2 {L) \\v £ \\ L 2 {L) ) , 
According to the trace Theorem and Sobolev's imbedding Theorem, we have 

Kilmer) < c \\ v e\\m/i{T) ^ c he\\m(p.), 
\\Ve\\l2(L) < c\\v E \\m/2(L) < c \\ v e\\m(n)- 

Hence, according to (l56l >. we get 

(57) 11^11^(0) < c\\^eU 2 \\ L 2(Y) + c\\l + d n U 2 \\ L 2(L)- 

Now we prove that ||V ; £' u 2||L 2 (r) — >• as e — ^ 0. We may assume that the system of cartesian coordinates 
(O, xi, x 2 ) is such that the origin O is one of the points A\ or A 2 and that T is locally above the xi-axis; see 
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Figure 4: T is locally the graph of a convex function, with a tangent to the X2-axis. 



Figure |4] Since Q is convex, there exist 5 > and two constants a > and /3 such that for all x\ G (0, 5), 
T is the graph of a convex function / of x\. For our choice of ip £ , since supp ip £ = [0, /3 e ], we have the 
estimate 



According to lfl2l |2D and our previous remarks in section |4~T1 we have U2 = y/r cos(#/2) + Uqo, with 
Uoo G H 2 (Q), and (r, ^) are the polar coordinates defined previously with origin 0. Thus there exists a 
constant c such that 

1 1*2 1 < c-y/r cos(#/2) 

in a neighborhood of with G (0,7r/2). Indeed, Uqo is # 2 therefore C 1 in a neighborhood of and 
then has an expansion of the form: Uqo = c s r + o(r), as r — > 0. Note that r = \J x\ + x\ and thus 
r = \J x\ + f(xi) 2 on T. Then 



The function / is convex and /(0) = 0, thus /' > for e small enough. Since the boundary T is tangent to 
the (0x2 ) axis, we have 





f'(x\) — > oo as xi — ^ 
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Thus, for e > small enough 



5«2||i2(r) < c^ e (0) ^ f(xi)f'(xi)dxi 



< 



opM (/(/3 e ) 2 ) 1/2 = c^(0)/(/3 e ). 
Since /(xi) — > as xi — > 0, we may choose ^> £ (0) and @ e in order to obtain ip £ (0)f(f3 £ ) ~ > as e — > and 

(58) ||^e«2[Ua(r) -> as e 0. 

Then, in view of (I57T ). we may deduce that ||u e ||ij-i(n) is bounded for the appropriate choice of ip £ . Conse- 
quently, ||u e [|.z, 2 (r) an d [l v e||z 2 (£) WQ a l so bounded. Using d55l ), we may also write 

&(0)KHi a(i) = / (Ve)V 6 < / K)Ve 
J L J dfl 

(59) < ||^ e W2||L2(r)lke||i 2 (r) + II 1 + d n U 2 \\ L 2^\\v £ \\ L 2^ L y 

Since ^> e (0) — )• oo as e — >• and all terms in d59l > are bounded, we necessarily have 

IMIl 2 (l) as e 0. 
Finally going back to (1561 and using the previous results, we obtain 

IKII-ff 1 ^) Oase 0, 

and this proves U2, £ — > U2 as e — > 0, in □ 

The following theorem gives a mathematical justification of the numerical scheme implemented in sec- 
tion[6]to find the solution of the free Bernoulli problem (J 7 ), based on the use of a penalized functional J £ 
defined by 

(60) J e (fi) = [ (u 2i£ - Ul ) 2 , 



where u\ is the solution of (I30l)-(f32l) and u 2j£ is the solution of (l47T)-(|49l). 
Theorem 3. One has 

lim inf(J £ (0)- J(O))=0. 
e^0 QeO 

Proof. The main ingredient of this proof is the result stated in Proposition [2] Indeed, this proposition yields 
in particular the convergence of U2, £ to 112 in L 2 (Q), when Q is a fixed element of O. It follows immediately 
that 

J e (0) -> J(n), ase^ 0. 
Let us denote by the solution of the free Bernoulli problem (J 7 ). Then, we obviously have 



mf jju) < Jjn*). 



Then, going to the limit as e — > yields 



< lim inf JJQ) < lim JJQ*) = J(Q*) = 0. 

e->0 fleO e^0 



□ 
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Remark 2. Theorem \3\ does not imply the existence of solutions for the problem inf{J e (0), Vt E O} and 
the following questions remain open: ( i) existence of a minimizer 0* for this problem, ( ii) compactness of 
(Vt*)for an appropriate topology of domains. These problems appear difficult since to solve it, we probably 
need to establish a Sverak-like theorem for the Laplacian with Robin boundary conditions and some counter 
examples (see e.g. suggest that this is in general not true. 

Nevertheless, if ( i) and ( ii) are true, Theorem \3\ implies the convergence as e — > 0, of f2* to £1*, the 
solution of CD)-©- 

5 Shape derivative for the penalized Bernoulli problem 

In order to stay in the class of domains O, the speed V should satisfy 

(61) V(x) = Vx G K, 

(62) V(x) ■ n(x) < Vx G L. 

Condition (I6TT) will be taken into account in the algorithm, and (1621 will be guaranteed by our optimization 
algorithm. We have the following result for the shape derivative dJ £ (£l; V) of J £ (£l) 

Theorem 4. The shape derivative dJ £ (£l; V) of J £ at 0, in the direction V is given by 

dJ £ (O; V) = J (Vpi • Vui + Vp 2 • Vn 2 , e + p 2 n + (ui - u 2y£ f) V -ndT, 

+ J (Vpi • Vui - Vp 2 ■ Vu 2 , £ ) V -ndL, 

where H is the mean curvature ofT and pi, p 2 are given by (172l) - (173l) and (174l) - (176l ). respectively. 
Proof. According to [6l[T5]|27]], the shape derivative of J £ is given by 



(63) dJ e (n;V)= [ 2(«!-« 2 ) «-«')+ f (m -u 2 . £ ?V ■ n, 

Jn Jan 

where u[ and u' 2 £ are the so-called shape derivatives of u\ and u 2 , respectively, and solve 

(64) -Au[ =0 in Q, 

(65) u' x = on K, 

(66) u'i = — d n u\V ■ n ondQ\K, 

(67) -Aua )£ = in n, 

(68) u' 2 e = on K, 

(69) e = — d n v,2, e V • n on L, 

9 n n 2 e + ip £ u' 2jE = div r (V ■ nV r n 2>e ) 

(70) -HV-n- i) £ d n u 2:£ V ■ n on T, 



where H denotes the mean curvature of T, and Vr is the tangential gradient on T defined by 

Vpu = Vu — (d n u)n. 
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Note that u' x and u' 2 £ both vanish on K, indeed, K is fixed due to (I6TT) which follows from the definition of 
our problem and of the class O. Further we will also need 

(71) d n u' 2 E = divr(V ■ nV r u 2 ,e) - UV ■ n on T, 

which is obtained in the same way as ( TTOl ). We introduce the adjoint states p\ and p 2 



(72) 


-Api = 2(m - u 2 , e ) 


in 17 


(73) 


pi = on Oft, 




(74) 


-Ap 2 = 2(m - U 2 , e ) 


in SI 


(75) 


p 2 = on L U K, 




(76) 


d n P2 = on T. 





Note that pi and p 2 actually depend on s although this is not apparent in the notation for the sake of read- 
ability. Using the adjoint states, we are able to compute 

/ 2(m - U2, e )u'x = / -Apiu'i 

Jq Jn 

= / -Au[pi - / d n piu[ - pid n ui 
Jn Jan 

d n piu[ 

an\K 

d n p\d n u\V ■ n. 

ldU\K 

Observing that Vpi = d n p\n and Vu\ = d n u\n on dVl \ K due to (1321 and (1731 we obtain 



(77) / 2(m - u 2)£ )u[dx = / Vpi • VuiF • n. 



For the other domain integral in (1631 we get 

/ 2(m - u 2:£ )u' 2 e = / -Ap 2 u' 2jE 
Jn Jn 

= / -^ 2z V2- \ {d n p 2 u' 2>£ -p 2 d n u' 2e ). 
Jn Jan 

At this point we make use of (l67"T)-(r7TI) and we get 

/ 2(m - U2,e)u' 2 £ = / p 2 (divr(V • nVrU2,e) - UV ■ n)dT + / d n p 2 d n u 2i£ V ■ n dL. 
Jn Jr Jl 

Applying classical tangential calculus to the above equation (see [27, Proposition 2.57] for instance) we 

have 

/ 2(«i - U2,e)u 2 £ = ~ (VrP2 ■ V r ti2,e^ • n - P2HV ■ n)dT + I d n p2d n u 2 , £ V ■ n dL 
Jn Jt Jl 

= -J (Vp 2 ■ Vu 2)£ V ■ n - p 2 nV ■ n)dT + J Vp 2 - Vu 2 , £ V ■ n dL, 
and the proof is complete. □ 
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6 Numerical scheme 



6.1 Parameterization versus level set method 

For the numerical realization of shape optimization problems, the main issue is the representation of the 
moving shape Q. Several different techniques are available: for our purpose, the most appropriate methods 
would be parameterization and the level set method. In the parameterization method for two-dimensional 
problems, curves are typically represented as splines given by control points = (£1^, £2,fc)> k = 0, .., m 
with m G IN*. The coordinates of these control points then become the shape design variables. In the level 
set method, the boundary of the domain in H N is implicitely given by the zero level set of a function in 
R iV+1 . Parameterization methods are the easiest to implement if the topology of the domain £1 does not 
change in the course of iterations, whereas the level set method is more technical to implement but thanks to 
the implicit representation, it allows to handle easily topological changes of the domain, such as the creation 
of holes or the merging of two connected components. 

For instance, in |@l|22l, the level set method is used to solve Bernoulli free boundary problem where the 
number of connected components is not known beforehand. In our case, we are solving the free boundary 
problem (J 7 ) in the class O of convex domains, thus the domains only have one connected component and 
the topology is known. In this case it is better to opt for the parameterization method which is easier to 
implement and lighter in terms of computations. 

The free boundary T C dQ is represented with the help of a Bezier curve of degree m G IN*. Let 



x(s) = (xi(s),x 2 {s)) 



s G [0, 1] 



be a parametric representation of the open curve T and let 



6c = (£l,fc,£2,k)j k 



0, .., m 



be a set of m + 1 control points such that the parameterization of T satisfies 



m 



(78) 




where 



(79) 




and (™) are the binomial coefficients. The geometric features such as the unit tangent t(s), unit normal 
n(s) and curvature %{s) are easily obtained from the representation (|78T >. Indeed we have 



(80) 



r(s)=x'(s)/\x'(s)\ 



with 



in 



(81) 




k=0 
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The coefficients B' k (s) are derived from d79l 



(82) B' 



k,m 



S 



(T) h fc ^ (1 " S ^~ H ^ + ( k ~ ~ s ) m " fc-:l l{ft<m-l} 



Since n(s) • r(s) = 0, we deduce the expression for the unit normal n(s) 

Y,k=0 Blm( S )£k 



(83) n{s) 

J2k=0 ^'k,m( S )^k 

with & := (S,2,k, —£,i,k)- The curvature %(s) is obtained with the help of formula 

(84) t'{s) =%{s)n{s). 
Thus we take 

(85) n(s) = t'(s) ■ n{s). 
Remark 3. According to (l80l >, (|8TT > arad (l82l , we obtain 



(86) r(0)= ,^ | T , r(l) ^ * ro " 1 



ICl - Col |£m _ Cm-ll 

Thus, in order to create a curve which is tangent to the axis {x\ = 0}, we need to take ^q, £i and (,m-i,£,m 
on {x\ = 0}. 

6.2 Algorithm 

For the numerical algorithm we use a gradient projection method in order to deal with the geometric con- 
straint C IR+; see the textbooks |[20ll25l for details on the method. A solution for dealing with the shape 
optimization problems with a convexity constraint is to parameterize the boundary using a support function 
w. If one uses a polar coordinates representation (r, 6) for the domains, namely 

fl w := \ (r, 9) G [0, oo) x E; r < 



w{6) 

where w is a positive and 2n -periodic function, then Q w is convex if and only if w" + w > 0; see ||2"3l for 
details. However, in our case, the convexity constraint for Q, is not implemented (i.e. we relax this constraint) 
for the sake of simplicity, but the convexity property is observed at every iteration and in particular for the 
optimal domain if the initial domain is convex. Moreover, Theorem 6.6.2 of [ 15'] may be easily generalized 
in our case and guarantees the convexity of the solution of the free boundary problem (J 7 ) even if the 
convexity hypothesis were not contained in the set O. 

We will denote by a superscript (I) an object at iteration /. The algorithm is as follows: we are looking 
for an update of the design variable of the type 

(8V) +ad^), 

where P stands for the projection on the set of constraints and a is the steplength which has to be determined 
by an appropriate linesearch. In our case, the constraint is Q C R^, which implies the constraint 

(88) ziO)>0, VsG[0, 1]. 
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In view of (|78T ). it is difficult to directly interpret the constraint (1881) for individual control points We 
choose therefore to impose the stronger constraint 



£i, fc >0, Vfee{l,..,m}. 



(89) 



for the control points. Constraint d89l is stronger than (I88I ). indeed, on one hand there might exist a such 
that ^i fc < while (I88T ) is still satisfied, but on the other hand, condition (l89l implies (188T ). However, in 
our case, the tips x(0) and x(l) of T are moving and the constraint should not be active for the points of 
T on the optimal domain. With (|89l we only guarantee that the domain stays feasible, i.e. f2 G IR+ for all 
iterates. In view of Remark[3l we also impose 



£2,0 — £2,1 — £2,m-i — 6 



2,777 







Al+i) 



max U^ k + ad^ k ,0 



in order to preserve the tangent to the axis {x\ = 0} at the tips of V. Therefore, for A; = 0, .., m, is 
updated using, 

(90) 
(91) 
(92) 
(93) 



?2,fc — ?2,fc ^ ««? 2 ,fc' 



"S2,m-1 — a ?2,m 



0. 



The link between the perturbation field V and the step d£f, is directly established using 

m 

(94) V(x(s)) = J2 B k,m( s ) d tk- 

k=0 

Thus, with a shape derivative given by 

(95) dJ E (Q;V)= [ VJ e (x)V(x)-n(x)dT(x) 

Jdn 



as in Theorem 01 we obtain using (l94l and ((95 



dJ £ (ft; V) = VJ e (x(s))V(x(s)) ■ n(s)\x'(s)\ ds 



VJ £ (x(s)) 



,fc=0 



n(s)|x'(s)| <is 



Vd^fc- / VJ B (z(s))5 fc>m (s)n(s)|x , (s)|ds. 



Thus, a descent direction for the algorithm is given by 



(96) 



VJ £ (x(s))B km (s)n(s)\x'(s)\ ds, 



and we obtain 
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and the update is then performed according to d90l)-(l93l). The step a is determined by a line search in the 
spirit of the gradient projection algorithm ll20l : a step is validated if we observe a sufficient decrease of the 
shape functional J £ measured by 

m 

Je{^)-u^)<-lY,\€ +l) -^\\ 

k=\ 

where | • | denotes the Euclidian distance. The line search consists in finding the smallest integer a (the 
smallest possible being a = 0) such that 

a = nr] a , 

where n and 7] < 1 are user-defined parameters. To stop the algorithm, we use the following stopping 
criterion: we stop when 

_Ai) I < _ _ t(o)i 

where r r is a user-defined parameter. 

7 Numerical results 

For the numerical resolution we take m = 40 control points We discretize the interval [0, 1] for the 
parameterization x(s) using 400 points. The domain K is chosen as 

K = {0} x [0.5-ki,0.5 + ki], 

with ki rj 0.129. The initial domain L is chosen as 

L = {0} x [0.5 - k 2 , 0.5 - m] U [0.5 + «i, 0.5 + k 2 ], 

with K2 ~ 0.233. We use the Matlab PDE toolbox to produce a grid in 17 and solve ui,U2 t£ >Pi,P2 using 
finite elements. The geometric quantities such as tangent, normal and curvature are computed using (l80l >- 
(|8TT ). (|83T > and (|85T >. respectively. We initialize the points ^ by placing them evenly on a half-circle of center 
{0} x {0.5} and radius 0.3, except for the two first £q, £i and two last points £ m _i, £ TO which have to lay on 
the axis {x\ = 0} as mentionned earlier. We choose fi = 10, rj = 0.5 for the line search and r r = 5 x 10~ 4 
for the stopping criterion. For the penalization we use (l46l) and choose e = 10 _1 and q = 4. 

The algorithm terminated after 220 iterations. The results are given in Figures [5] to |7] In Figure [51 the 
two states u\ and U2, £ as we U as tne two adjoint states pi and p 2 are plotted. The difference between u\ and 
u 2)£ in tne fi na l domain £lfi na i is plotted in Figure [6l along with the residual J £ (Q) given by (l60l) . In Figure 
|7J the initial and final boundaries are plotted in blue and red, respectively, while the set of control points of 
the curve T is plotted in green. We observe that the optimal domain is symmetric as expected from section 
13.11 The optimal set L j ina i is given by 

Lfinal = {0} X [0.5 -K /i5 0.5-Ki]U [0.5 + «i, 0.5 + K/j]. 

with Kji Rj 0.2342. The value of J £ on the initial domain is 

Jei^initial) ~ 2.6 X 10" 3 , 
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Figure 5: Solutions u\ (top left), U2, £ (top right), p\ (bottom left), p2 (bottom right) in the optimal domain, 
and the value of J £ on the final domain is 

Jei^final) ~ 3.3 X 1(T 8 , 

as may be seen in Figure [6l Therefore, the shape functional J £ has been significantely decreased and is close 
to its global optimum. 
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Figure 6: Difference u\ — U2 in the optimal domain (left), residual J £ (right). 
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Figure 7: Final boundary T (red), initial boundary T (blue), control points (green). 
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